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ABSTRACT 

This  report  presents  a  multivariable  controller  design  for  a  generic  missile  travelling  at 
hypersonic  speed.  The  controller  is  based  on  the  theory  of  optimal  control,  and  the  solution  is 
provided  by  a  feed  back  controller  known  as  a  Li  near  Quadratic  Regulator  (LQR)  for  tracking  the 
demands  (lateral  acceleration  and  roll  rate)  as  closely  as  possible  while  keeping  the  actuator 
efforts  (for  deflecting  ai  leron,  elevator  and  rudder)  smal  I .  I  n  the  computer  simulation,  the  missile 
model  is  subject  to  various  hypersonic  flight  conditions  and  gain  scheduling  is  exercised  as  a 
function  of  dynamic  pressure.  The  inti  al  study  assumes  perfect  knowledgeof  all  the  states  for  the 
controller  design.  Subsequently,  a  state  observer  (kal  man  filter)  is  introduced  into  the  design  to 
reflect  a  more  practical  synthesis  where  some  of  the  states  need  to  be  estimated. 
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Design  of  Optimal  Flight  Controller  for  Generic 
Linearised  Missile  Model  in  Hypersonic  Regime 


Executi  ve  Su  mmary 


The  design  of  a  missile  autopilot  controller  that  provides  a  fast  and  stable  response  to 
guidance  commands  and  also  maintains  the  performance  throughout  the  entire  flight 
envelope,  poses  a  challenge  to  the  flight  control  engineer.  This  report  considers  the 
autop  i  I  ot  d  esi  gn  for  a  generi  c  cru  ci  form  mi  ssi  I  e  uti  I  i  si  ng  the  generi  c  6D  oF  ai  rf  rame  mod  el 
developed  in-house[Faruqi  2007]  and  the  linear  quadratic  regulator  (LQR)  technique. 

The  control  designs  are  presented  for  tracking  the  desired  lateral  accelerations  and  roll 
rate,  where  the  control  inputs  are  aileron,  elevator  and  rudder  deflections.  The  optimal 
LQR  design  aims  to  find  the  feedback  gain  that  minimises  a  particular  cost  function.  In 
this  case  it  is  the  weighted  sum  of  energies  associated  with  tracking  error  and  control 
si  gnal .  N  ote  that  these  two  cost  terms  areconfl  i  cti  ng  (i  ,e„  faster  tracki  ng  I  eads  to  i  ncreased 
strain  on  the  actuator)  and  the  role  of  the  control  gain  is  simply  a  trade-off  between  the 
response  time  and  control  effort. 

Wesubject  thevehideto  two  hypersonic  flight  conditions  and  show  that  a  gain  scheduling 
may  be  performed  using  thedynamicpressureasthesoleschedulinginputThecontroller 
design  extends  to  output  feedback  control  where  a  state  observer  is  included  in  the 
feedback  loop. 
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1.  Introduction 

The  design  of  a  missile  autopilot  controller  which  provides  a  fast  and  stable  response  to 
guidance  commands,  and  which  maintains  the  performance  throughout  the  entire  flight 
envel  ope,  poses  a  chal  I  enge  to  the  fl  i  ght  control  engi  neer.  Thi  s  report  consi  ders  the  autopi  I  ot 
design  for  a  generic  cruciform  missile  utilising  the  generic  6D0F  airframe  model  developed 
in-house  [Faruqi  2007]  and  the  I  inear  quadratic  regulator  (LQR)  technique. 

Classical  control  techniques  usually  involve  loop  shaping  approaches  -  frequency  domain 
techniques  such  as  bode  and  Nyquist  plots.  Single-Input  Single-Output  (SI SO)  systems  are 
the  primary  focus  of  classical  controls.  M  ultiple-l  nput  M  ultiple-Output  (M I  MO)  systems  can 
be  controlled  through  classical  control  techniques,  but  to  do  so  usually  requires  that  each 
input  and  output  can  bedecoupled.  In  this  case  the  system  can  be  treated  asabunch  of  SI  SO 
systems  and  each  controller  can  bedesigned  without  consideration  of  the  others. 

Controllers  are  filters  that  shape  the  feed  back  error  to  achieve  the  desired  system  response. 
The  input  to  a  controller  is  the  error  between  the  desired  and  measured  system  outputs.  The 
controller's  output  is  a  filtered  version  of  that  error.  Controller  design  is  a  matter  of  selecting 
poles,  zeros,  and  gains.  Poles  provide  integral  action  which  is  slow  responding  to  errors  but 
drives  the  step  response  to  a  zero  steady-state  error.  Zeros  are  responsible  for  derivative 
action  which  is  fast  but  can  amplify  noise  and  leads  to  a  non-zero  steady-state  error.  Gains 
provide  instantaneous  response  to  errors. 

When  applied  to  the  flight  control,  the  classical  approach  leads  to  a  three-loop  autopilot 
[Neslineand  Nabbefeld  1984],  where  longitudinal,  lateral  and  roll  dynamics  are  considered 
independently  as  mentioned  above.  This  approach  is  valid  for  a  limited  range  of  flight 
conditions.  Furthermore,  loop  shaping  can  bevery  difficult  for  MIMO  systems. 

The  state-based  modern  control  design  allows  for  a  more  convenient  and  compact  way  to 
model  and  analyseMIMO  systems,  where  the  inter-axial  coupling  effects  are  included.  This 
approach  can  be  applied  to  more  extreme  flight  regimes,  non-symmetrical  airframes  with 
non-symmetrical  mass  distribution,  and  adaptive  control. 

The  optimal  LQR  design  aims  to  minimise  or  maximisea  particular  cost  function.  In  our  case 
it  is  the  minimisation  of  the  combined  energy  of  tracking  error  (response  speed)  and  control 
si  gnal .  N  ote  that  these  two  cost  terms  are  confl  i  cti  ng  (i  .e.  decreasi  ng  one  i  ncreases  the  other) 
and  the  role  of  control  gain  is  simply  the  trade-off  bdween  the  response  time  and  control 
effort. 

We  initially  assume  that  all  24  states  are  precisely  known,  and  design  a  full  state  feedback 
controller.  We  then  present  the  stability  issue,  steady  state  gain,  effect  of  poles  and  zeros  on 
step  response,  and  gain  scheduling  for  the  two  flight  envelopes  associated  with  the  hypersonic 
air-to-ground  attack  scenarios.  We  then  move  on  to  the  observed  feedback  control  design, 
where  only  6  inertial  measurement  unit  (IMU)  output  variables  are  available,  and  a  Kalman 
f  i  I  ter  i  s  used  to  esti  mate  the  24  state  vari  abl  es  from  the  6  measu  rement  vari  abl  es.  Repl  aci  ng  the 
true  state  variables  with  the  estimates,  we  repeat  some  of  the  previous  computer  simulations 
and  i  nvesti  gate  the  effect  of  the  state  observer  on  the  cl  osed  I  oop  response. 
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2.  State  Feedback  Control  System 

In  this  section,  we  present  the  closed-loop  optimal  control  of  linear  systems  with  quadratic 
performance  index.  This  leads  to  the  linear  quadratic  regulator  (LQR)  system  dealing  with 
state  regulation  (keeping  the  state  around  zero).  If  the  requirement  is  altered  to  make  the 
state  follow  or  track  a  desired  trajectory,  then  the  problem  becomes  that  of  linear  quadratic 
tracking  (LQT)  or  alternatively  the  servomechanism  problem. 

2.1  Formulation 

Thelinearised  dynamical  system  is  governed  by  the  foil  owing  differential  equations 

x(t)  =  F(?)x(r)  +  G(r)u(?)  +  w  (t)  (1) 

and  the  measurement  equations  are  given  as 

z(t)  =  J(t)x(t)  +  v(t).  (2) 

where 

x  i  s  a24x  l  system  state  vector, 

uisa  3xl  (servo)  input  vector  =[£,?/,  £]r , 

z  is  a  6x  l  measurement  (accelerations  and  body  rates)vector  =  [ax,ay,az,p,q,r]T , 

w  i  s  a  24  x  l  process  noi  se  vector, 
v  i  s  a  6  x  l  measu  rement  noi  se  vector, 

F  i  s  a  24  x  24  system  coeffi  ci  ent  matri  x, 

G  i  s  a  24  x  3  control  coeffi  ci  ent  matri  x, 

J  is  a  6x24  measu  rement  coeffi  dent  matrix. 


The  state  vector  x  contai  ns  24  states: 


(u,v,w,p,q,r,ax,dx,a  d  az,dz,p0,p0,q0,q0,r0,r0,^0,i0,Ti0,ii0X0,C0)  where  the  first  six 


elements  are  3  velocity  components  and  3  angular  rate  components,  the  next  12  elements  are 
the  accelerometer  and  gyroscope  outputs  and  their  rates,  and  the  remaining  6  elements  are 
the  control  surface  deflections  and  their  rates).  The  derivation  and  detailed  descriptions  of 
the  above  state  space  model  are  given  in  [Faruqi  2007].  Figure  1  shows  the  step  responses  of 
the  24  states  when  ay  demand  =  lms  2 . 
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time  (s)  time  (s)  time  (s)  time  (s)  time  (s)  time  (s) 


Figurel:  Step  response  of  24  states 

Our  objective  is  to  control  the  system  in  such  a  way  that  the  output  z (t)  tracks  the  desired 
output  ii(0  as  closely  as  possible  during  the  interval  [r^^  jsuch  that  the  cost  function 
shown  below  is  minimised.  Let  us  define  the  error  vector  as  e(t)  =  x\(t)-z(t) ,  and  choose  the 
performance  i  ndex  as 

C  =  -  er  (t  f  )S  (t  f  )e(f  f)  +  -\tf  [er  (0Q(0e(0  +  u  r  (0R(0«(01*  (3) 

2  2  J'o 

where 

S  isa  6x6symmetricpositivesemi-definiteweighting matrix, 

Q  isa  6x 6 symmetric positivesemi-definiteweighting  matrix, 

R  isa  3x3  symmetricpositivedefiniteinputweighting  matrix. 

From  now  on,  we  omit  t  for  brevity  except  for  the  initial  and  final  time. 

Taking  the  variational  approach  as  in  [Sage  1969;  Naidu  2002],  the  solution  to  the  optimal 
feedback  control  problem  is  given  as: 

u  =  [-R~'GrP]x  +  [R_1Gr]§  (4) 
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where 

P  =  -PF-FrP  +  PGR  ‘GrP- JrQJ  (Differential  Riccati  Equation)  (5) 

P(f/)  =  iT (t f)S(t f)i(t f)  and  (6) 

4  =  -[F-GR  GrP]r^  -  JrQr]  (7) 

=  (8) 

Given  the  final  conditions  (6)  and  (8),  the  differential  equations  (5)  and  (7)  can  be  solved  by 
integrating  backwards  in  time  (from  tf  to  t0 ).  Note  that  the  Riccati  equation  (5)  can  also  be 
solved  analytically  [Anderson  and  Moore  1989;  Naidu  2002], 

By  setting  tf  =  oo,  the  problem  becomes  the  infinite  time  (steady-state)  case,  and  the 
performance  i  ndex  si  mpl  ifies  to 

C  =  -  fV(0Q(0e(0  +  ur(0R(0u(0]*  (9) 

2  «o 

The  steady  stateP  isobtained  by  solving  theContinuous  Algebraic  Riccati  Equation  (CARE) 

0  =  -PF-FrP  +  PGR  GrP-JrQJ  (10) 

And  the  steady  state  %  is  solved  by  siting  \  =  0  and  is  given  by 

§  =  -[(F-GR  *GrP)r]  1  JrQn  (11) 

We  can  use  the  analytic  solution  [Anderson  and  Moore  1989;  Naidu  2002]  or  a  matlab 
toolbox  function  "care"  as  shown  below. 

P  =care(F,G,J7QJ,R) 

H  aving  found  the  stabilising  sol  utionsP  and  % ,  the  optimal  feedback  gain  Kc  is  given  as 
Kt  =R  GrP  (12) 

Adding  the  input  gain  term,  the  solution  to  the  control  input  is  written  as 

u  =  -Kfx  +  K^  (13) 

where  Kj  =R~‘Gr  (input  gain  matrix)  (14) 
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Weexpress  c  as  (refer  to  Figure  2), 

(15) 

where 

A  =  -(F-GR  GrP)r  =  -(F-GKf)r  (16) 

B  =  JrQ  (17) 

Note  that  for  the  linear  quadratic  regulator  (LQR)  problem,  the  input  gain  term  (11)  is 
omitted  and  J  term  in  (10)  also  disappears.  The  transient  behaviour  and  stability  of  the 
demand  tracking  problem  is  the  same  as  the  regulator  problem,  therefore  the  stability 
analysis  is  performed  in  LQR  context. 

2.2  Stability  of  LQR  system 

Ignoring  the  input  term,  the  dosed  loop  differential  equation  becomes, 

x(r)  =  (F-GKc)x  (18) 

In  order  for  F -GKf to  be  a  stable  matrix,  the  second  method  of  Lyapunov  [Ogata  1995] 
states  that  there  must  exi  st  a  positive  defi  nite  matrix  P  that  sati  sfi  es  thefol  I  owi  ng  equati  on 

[(F-GKe)rP  +  P(F  -GKC)]  +  X  =  0  (19) 

whereX  is  some  positive  definite  symmetric  matrix. 

We  know  Q  is  positive  semi-definite  and  R  is  positive  definite,  and  they  are  both  diagonal, 
hence  w e  may  say  Q  +  K Tc  RK f  i  s  posi ti ve  d efi  n i te.  We  i  ntrod u ce  a  matri x X  and  set  i t  as 

X  =  Q  +  K^RK  (20) 

Substituti  ng  X  =  Q  +  K^RKf  i  nto  (19),  we  get 

(F  -  GKf  )rP  +  P(F -GKc)  +  (Q  +  K^RKc)  -  0  .  (21) 

Expanding  the  above  equati  on  leads  to 

FrP  +  PF  - K^GrP  - PGKC  +Q  +  K^RKC  =  0  .  (22) 

The  optimal  feedback  gain  according  to  the  LQR  design  isKc=R_1GrP  as  in  (12). 
Substituting  R  GrP  for  Kf  in  (22)  leads  to  thefollowing  stability  condition 

FrP  +  PF  -  PGR  GrP  +  Q  =  0  ,  (23) 
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This  is  the  algebraic  Riccati  equation  in  the  LQR  problem  that  was  used  earlier  to  compute 
P  which  was  utilised  to  generate  Kc  and  u  .  Since  the  necessary  stability  condition  (i.e.,  the 
second  method  of  Lyapunov  [Ogata  1995])  leads  to  the  Riccati  equation,  we  may  say  that  the 
proposed  closed-loop  system  is  asymptotically  stable. 


Figure  2:  LQT  design  with  G  aussian  error  terms  and  state  estimator  (output  feedback) 


Figure  3:  E  rror  free  LQT  desi  gn  with  full  state  feedback 

The  proposed  linear  quadratic  tracker  (LQT)  design  with  the  feedback  and  input  gains 
Kt  and  ,  and  state  observer  is  shown  in  Figure  2.  However,  the  majority  of  the 
simulations  in  this  report  are  based  on  the  full -state  feedback  scheme  as  in  Figure  3. 

2.3  Steady  State  Gain 

Substituting  the  optimal  input  and  feedback  gain  terms  into  (1),  we  get  the  foil  owing  closed- 
loop  differential  equation 


x  =  [F  -GKc]x  +  [GKj((GKc  -F)r)  1  JrQ]ii . 


(24) 
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A  t  steady  state,  x  =  0,  and  the  analytical  solution  for  steady  state  x  isgivenas 

x  =  [GKC  -F]-‘GK1[(GKC  -F)r]  1  JrQrj  .  (25) 

By  setting  R  =  I3x3 ,  we  get  =  Gr  (see  (14)),  and  Equation  (25)  becomes 

X  =  [GKC  -F]  GGr[(GKc  -F)r]  1  JrQii  (26) 

The  measurement  equation  is 


z  =  Jx 

-  J(GKC  -F)  GGr((GKc  -F)r)  JrQii 
=  [J(GKC  -F)-1G][J(GKC  -F)-‘G]rQii 

=  Tssil  (27) 

The  steady  state  gain  matrix  Tss  =  [J(GKC -F)_1G][J(GKe -F)  *G]rQ  is  symmetric  and 
maps  the  desired  output  (or  demand)  r\  to  the  actual  outputz  in  steady  state. 


Typically  Tss  appears  as  below 

"0  0  0  0  0  O' 

0  1  0  0  0  0 

0  0  1  0  0  0 

T  « 

ss  0  0  0  1  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 


The  non-zero  diagonal  elements  of  Tss  represent  steady  state  gains,  ay/ay  demand  , 

«z/«z  demand  and  P0  / ^demand  in  descending  order.  The  remaining  terms  are  either  zeros  or 

close  to  zero.  The  non-zero  diagonal  elements  approximate  to  1  only  if  Q  is  large.  It  was 
observed  that  as  Q  gets  smaller,  the  diagonal  terms  get  smaller  as  well.  In  order  to  explain 
this  observation,  wepresent  a  simplescalar  differential  equation. 


x  =  fx  +  u  (28) 

z  =  x  (29) 

Since  g  =  r  =  1 ,  the  gains  become  Kc-p  and  Kl  =  1 ,  then  the  scalar  version  of  (25)  becomes 

*  =  (f-p)x  +  — -  (30) 

p-f 
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At  steady  state  (ie„  x  =  0 ), 


qrj 


(. p-fY 


(31) 


w  here  q  i  s  a  seal  ar  wei  ghti  ng  coeffi  ci  ent. 
Then  the  measurement  output  becomes 


Z  =  X  —  ■ 


q> i 


( p-f ) 


2  ' 


(32) 


Expanding  the  denominator,  the  gain  can  be  shown  as 


q 


q  p--2pf  +  f~ 


(33) 


Now  solving  the  scalar  Riccati  equation,  2pf  - p2  +q  =  0  ,  we  obtain  the  optimal  q  as 
q  =  p2  -2pf ,  therefore  the  above  gain  can  be  expressed  as 


z  _  q 
rj  q  +  f2  ' 


(34) 


It  is  clear  that  as  q  »  f2 ,  the  gain  —  ->  1 ,  and  as  q  «  f2 ,  —  ->  «  0  .  So  one  may  say  in 

q  q  f 

general  that  low  Q/R ratio  leads  to  the  steady  output  gain  that  is  smaller  than  one.  To 
ensure  the  unity  static  gain,  the  demand  needs  to  be  scaled  before  entering  the  system.  The 
lateral  acceleration  and  roll  rate  demands,  r|(2) , tj(3)  and  n(4),  need  to  be  divided  by 
Tss (2,2)  Tss  (3,3)  and  Tss (4, 4) ,  respectively  before  entering  the  system. 

In  the  following  examples,  we  look  at  some  of  the  selected  runs  with  different  dynamic 
pressure  values,  and  examine  their  responses. 

2.4  Step  Responses  at  D  ifferent  0  perati ng  Poi  nts 

The  simulations  presented  here  are  based  on  two  flight  scenarios  -  high  and  medium  height 
launch  modes.  Given  the  altitude  and  Mach  number,  dynamic  pressure  q  is  computed  as 

q=\pV2  (35) 

where  p  and  V  are  air  density  and  speed.  The  dynamic  pressure  values  are  computed  at 
various  operating  points  along  the  flight  envelopes,  and  are  listed  in  Table  land  Table  2. 
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Table  1:  H  igh  Launch  Flight  Envelope 


Altitude(m) 

M  ach  N  o 

Velocity  (m/  s) 

Dynamic  Pressure 

30000 

0.8 

241 

521 

30000 

2 

604 

3256 

30000 

4 

1207 

13026 

30000 

6 

1811 

29308 

27500 

8 

2401 

76102 

25000 

8.2 

2448 

117281 

22500 

8.4 

2493 

181313 

20000 

8.6 

2538 

283452 

17500 

8.8 

2597 

440208 

15000 

9 

2656 

682947 

12500 

9.2 

2715 

1058492 

10000 

9.4 

2815 

1637048 

7500 

9.6 

2978 

2470386 

5000 

9.8 

3142 

3635414 

2500 

10 

3306 

5232914 

0 

10.2 

3451 

7386198 

Table 2:  M  edium  Launch  Flight  Envelope 


Altitude(m) 

M  ach  N  o 

Velocity  (m/  s) 

Dynamic  Pressure 

15000 

0.8 

236 

5396 

20000 

2 

590 

15330 

25000 

4 

1194 

27908 

30000 

6 

1810 

29308 

30000 

8 

2415 

52103 

30000 

10 

3018 

81411 

27500 

10 

3002 

118910 

25000 

10 

2985 

174422 

22500 

10 

2968 

256962 

20000 

10 

2951 

383250 

17500 

10 

2951 

568450 

15000 

10 

2951 

843144 

12500 

10 

2951 

1250581 

10000 

10 

2995 

1852703 

7500 

10 

3102 

2680540 

5000 

10 

3206 

3785313 

2500 

10 

3306 

5232914 

0 

10 

3403 

7099383 

The  open-  and  closed-loop  responses  at  the  selected  operating  points  are  simulated  and  the 
corresponding  responses  are  examined.  Firstly  we  present  a  quick  look  comparison  between 
the  open  loop  and  closed  loop  responses,  for  a  selected  case  where  the  altitude  is  30  km  and 
the  Mach  number  is  2. 
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ay 


Figure4:  Open  and  closed  loop  responses.  Open  loop  response  is  shown  in  blue 


In  Figured  the  open  loop  response  (blue)  is  a  lightly  damped  oscillation,  with  a  frequency  of 
2.29  rads"1.  This  corresponds  to  the  third  dominant  mode  shown  below.  All  of  the  poles  have 
negative  real  parts,  indicating  that  the  system  is  open-loop  stable. 


Eigenvalue 

Damping 

Freq.  (rad/s) 

-3.65e-005 

-6.26e-001 

-7.90e-002  ±  2.29e+000i 

1.00e+000 

1.00e+000 

3.45e-002 

3.65e-005 

6.26e-001 

2.29e+000 

Regarding  the  closed  loop  responses,  all  of  them  converge  correctly  to  the  demand.  Note 
that  the  increase  in  Q  results  in  a  faster  response.  From  (3),  it  is  apparent  that  a  higher  Q 
setting  puts  more  weight  on  reducing  the  sum  of  square  of  the  gaps  between^  and  z  curves 
in  Figured  hence  resulting  in  a  faster  response. 
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Step  Response 

n 

/  - — 

i 

> 


Figure5:  Larger  Q  means  moreweight  on  reducing  theintegral  squares  of  the  difference  between  tj 
and  z  curves.  T he  shaded  area  shrinks  as  Q  becomes  larger. 

We  now  look  at  four  examples  of  "high  altitude  low  speed",  "high  altitude  high  speed", 
"medium  altitude  high  speed",  "low  altitudehigh  speed"  cases. 

Firstly  for  the  high-altitude  low-speed  case,  the  step  response  is  sluggish  as  shown  in  Figure 
6,  due  to  low  dynamic  pressure.  It  takes  more  than  3  seconds  for  the  response  to  settle. 
Selecting  high  Q  tends  to  speed  up  the  convergence  at  the  expense  of  slightly  higher 
overshoot.  The  largest  overshoot  is  about  8  percent.  Note  that  it  is  important  to  make  sure 
that  the  servo  i  nput  val  ues  are  feasi  bl  e. 

According  to  Bryson's  rule  [Franklin  et.  al.  2002],  the  first  choice  for  the  weighting  matrices 
Q  and  R  is  to  select  diagonal  Q  and  R  with 

Q(i,i)  =  l/(maximum  acceptible  value  of  e(/))  where  i  =  {l _ 6}  ,  and 

R(y,  j)  =  1  /(maximum  acceptible  value  of  u (_/))  where  j  —  {l,  2, 3} . 

Bryson's  rule  is  often  just  the  starting  point  to  an  iterative  trial  and  error  design  procedure 
until  the  desired  closed  loop  performance  is  achieved.  In  this  report,  we  set  R  =  I  and  only 
vary  Qas  a  scalar  multiple  of  I.  The  optimum  Q  is  hand-picked  as  the  one  giving  the  fast 
responsewith  small  overshoot. 

It  is  interesting  to  note  that  there  is  a  sharp  peak  at  the  beginning  of  the  rise.  This  is  usually 
due  to  the  existence  of  dominant  closed-loop  zeros  in  the  left-half  of  thecomplexs-plane.  We 
deal  with  these  domi  nant  zeros  and  thei  r  effects  i  n  the  next  secti on. 

As  for  the  high-altitude  high-speed  run  (see  Figure  7)  the  response  is  faster,  reaching  the 
steady  state  within  0.3  seconds.  The  overshoot  is  increased  to  about  10%.  It  is  also  noticed 
that  the  performance  i  s  consi  stent  over  a  broad  range  of  Q  setti  ngs. 
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Another  interesting  observation  is  that  the  optimum  Q  value  is  smaller  when  the  dynamic 
pressure  is  higher.  The  increase  of  the  dynamic  pressure  causes  the  airframe  to  respond 
faster  without  havi  ng  to  i  ncreasethe  control  effort. 

Further  increase  in  Mach  number  together  with  altitude  reduction  (see  Figure  8)  results  in 
very  large  value  of  dynamic  pressure.  Again,  this  causes  the  optimum  Q  to  become  much 
smaller.  The  optimum  Q  is  chosen  as  0.0005x1  (green)  which  results  in  the  overshoot  of 
19%  and  the  settling  time  of  0.08  second.  Raising  Q  above  0.01  xi  introduces  oscillations. 

As  for  the  low-altitude  high-speed  run  (see  Figure  9),  the  dynamic  pressure  is  extremely 
high,  leading  to  a  very  small  Q.  Again  the  green  curve  is  picked  as  the  optimum  response 
(Q  =1  x  10~7I ).  The  settling  time  is  less  than  0.05  second,  and  the  overshoot  is  about  15%. 


ay 


Figure  6:  H  igh  altitude  and  low  speed  case  (Start  of  Boost).  30  km,  M  ach  0.8,  q  =  521 


0  0.05  0.1  0.15  0.2  0.25  0.3  0.35  0.4 

time  (s) 

Figure  7:  H  igh  altitude  and  high  speed  case  (End  of  Boost).  30  km,  M  ach  6,  q  =  29308 
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0  0.02  0.04  0.06  0.08  0.1  0.12  0.14  0.16 

time  (s) 

Figure8:  M  edium  altitude  and  high  speed  case  (D  ive).  15km,  M  ach  9,  q  =  682941 


0  0.01  0.02  0.03  0.04  0.05  0.06  0.07 

time  (s) 

Figure9:  Low  altitude  and  high  speed  case  (Impact).  Okm,  M  ach  10,  q  =  7386198 

2.5  Effect  of  Poles  and  Zeros  on  the  Step  Response 

We  select  the  high-altitude  low-speed  example  in  Section  2.4  to  take  a  closer  look  at  the 
initial  positive  bump  in  the  step  response.  The  transfer  function  ay  / ay  demand  is  numerically 

given  as 

_  1.99x1Q19s17  +---  +  2.94xlQ42 
fly  demand  S2‘  +  •  •  ■  +  2.94  X  1()42 

This  transfer  function  has  21  poles  and  17  zeros.  After  pole-zero  cancellations,  the  transfer 
function  is  truncated  to  have  only  four  most  dominant  poles  -2.66  +  2.24/ 
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and -135.9 +  146.06/,  and  two  most  dominant  zeros  0.0004  +  4.66/,  before  proceeding  to  the 
analysis  of  the  effect  of  dominant  zeros. 


Step  Response 


Figure  10:  Step  responses  of  truncated  and  original  transfer  functions:  Blue->  [4  poles], 
Green->  [4  poles,  2  zeros],  Red->[21  poles,  17  zeros]. 

Figure  10  indicates  that  the  truncated  model  with  4  dominant  poles  and  2  dominant  zeros  is 
a  good  approximation  -  the  green  and  red  curves  are  well  aligned.  It  is  clear  that  the  zeros 
are  the  most  dominant.  If  we  remove  the  zeros,  the  bump  disappears  (seethe  blue  curve  in 
Figure  10). 

We  want  to  examine  how  the  movement  of  the  zeros  along  the  real  and  imaginary  axes 
affects  the  initial  bump.  The  simulation  results  are  given  in  Figure  11  and  Figure  12. 

In  Figure  11,  as  the  zeros  move  away  from  the  origin  along  the  real  axis,  the  initial  bump 
becomes  less  prominent  by  means  of  the  trough  on  the  right  being  elevated,  and  eventually 
disappears  (The  cyan  curve  corresponds  to  the  case  when  the  zeros  become  less  dominant 
than  the  poles).  As  the  zeros  move  away  from  the  origin  along  the  imaginary  axis,  the  actual 
height  of  the  bumps  drops  and  eventually  disappears  (see  Figure  12). 

Revisiting  the  higher  speed  run  (Mach  6)  as  in  Figure  7,  the  initial  bump  is  significantly 
attenuated  when  compared  with  the  slow  speed  run  in  Figure  6.  For  Figure  7,  the  dominant 
poles  are  (-23.89  +  31. 14/, -127.24  ±129.25/)  and  the  dominant  zeros  are(-0.24±95.7/) .  These 
zeros  sit  much  further  from  the  origin  than  the  zeros  in  slow  speed  run  of  Figure  6,  resulting 
in  the  diminished  bump.  This  initial  bump  disappears  as  the  dynamic  pressure  further 
increases. 
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Step  Response 


Figurell:  Effect  of  moving  the  zero  along  the  real  axis.  The  real  part  of  zeros  are:  -0.004  (blue), 
-0.4  (green),  -1  (red),  and -4  (cyan). 
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Figure  12:  Effect  of  moving  the  zero  along  the  imaginary  axis.  The  imaginary  part  of  zeros  are: 
-  ±  2  (red),  ±4  (blue)and  ±8  (green). 

2.6  Why  Positive  Bump  when  Tail  Controlled? 

Itisgiven  in  [Garnell  and  East  1977]  that  the  I M  U  ismounted  0.5  m  fore  of  the  missile  centre 
of  gravity  (CG).  We  define  dx as  the  IMU  displacement  from  the  CG  along  the  missile's 
forward  axis(e.g„  dx  =0.5).  Considering  the  most  sluggish  response  in  the  high-altitude  low- 
speed  case,  they-axis  linear  acceleration  and  z-axis  (yaw)  angular  acceleration  are  plotted  in 
Figure  13.  The  initial  transient  portions  (red-circle)  are  zoomed  for  a  closer  look.  A  negative 
rudder  deflection  as  in  Figure  14  (ie.,  anticlockwise  when  viewed  from  above)  initially 
pushes  the  missile  to  the  left  before  the  clockwise  rotation  steers  it  to  right.  Therefore  the 
sensor  attheCG  should  experience  a  brief  negativey-acceleration  first  and  then  a  positive  y- 


Step  Response 

- imaginary  part  of  zeros  =-4 

—  imaginary  part  of  zeros  =-8 
- imaginary  part  of  zeros  =-2 
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acceleration.  This  phenomenon  is  known  as  non-minimum  phase  response  which  often 
ari  ses  among  the  tai  I  -control  I  ed  mi  ssi  I  es. 

The  positive  bump  in  Figure  13  is  the  result  of  the  accelerometer  being  placed  fore  of  the  CG, 
hence  sensing  the  lever  arm  effect.  The  positive  angular  acceleration  (right  plot  in  Figure  13) 
translates  to  the  positive  linear  acceleration  (e.g.,  ay  =  dj  )  that  happens  to  be  greater  than 
the  initial  negative  bump. 


ay  rdot 


Figure  13:  ajy  acceleration)  and  i-(yaw  acceleration)  responses.  The  plots  on  the  bottom  are  the 
enlarged  initial  responses. 


rudder  fin  deflection 


Figure  14:  Ru  dder  fi  n  defl  ecti  on 
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As  mentioned  earlier,  the  rudder  fin  deflection  is  mostly  negative  in  Figure  14.  The  fact  that 
the  negative  fin  deflection  leads  to  the  positive  yaw  motion  confirms  that  the  missile  is  tail- 
controlled.  Another  way  to  prove  that  the  missile  is  tail -control  led,  is  to  shift  the  IMU  back 
to  the  CG  (i.e.,  dr= 0)  and  show  that  the  initial  bump  becomes  negative.  In  Figure  15, 
without  the  lever  arm  effect,  the  initial  y-acceleration  bump  is  negative  (peaked  to  about  - 
0.5 m/s2)  as  expected.  Adding  this  -0.5 m/s2  to  lever-arm  effect  («0.8  m/s2)  gives  a  net 
acceleration  of  0.3  m/  s2.  This  value  is  similar  to  the  0.2  m/  s2  bump  in  Figure  13.  Furthermore, 
it  is  also  shown  that  shifting  the  IMU  totheCG  resultsin  one  of  thezeros  becoming  positive 
as  shown  Figure  16. 


ay  yaw  rate 


time  (s)  time  (s) 


ay  -  zoomed  yaw  rate  -  zoomed 


time  (s)  time  (s) 


Figure  15:  N  on-minimum  phase  response  for  y-acceleration  and  the  associated  angular  acceleration, 
when  thelM  U  is  at  theCG.  Bottom  row  is  thezooming  of  the  initial  bumps  in  thetop  row. 
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Figure  16:  Pole  and  zero  plot  when  the  IM  U  is  at  theCG.  Poles  are  shown  as  crosses  and  zeros  are 
shown  as  circles. 


2.7  Gain  Scheduling  -  Optimal  Weight  versus  Dynamic  Pressure 

Gain  scheduling  is  an  approach  that  uses  a  family  of  linear  controllers,  each  of  which 
provides  satisfactory  control  at  its  operating  point.  The  scheduling  variables  are  defined  as 
parameters  used  to  determine  the  operating  point  the  system  is  currently  at.  In  a  flight 
control  system,  they  are  usually  the  altitude  and  Mach  number  as  given  in  Tables  1-2.  This  is 
one  of  the  simplest  and  most  intuitive  forms  of  adaptive  control  (strictly  speaking  this  is  not 
adaptive  as  the  system  does  not  have  self-organising  features). 

For  a  given  dynamic  pressure,  we  iterate  the  step  response  simulations  over  various  Q 
settings  and  hand-pick  the  optimum  Q.  The  optimal  Q  values  are  then  plotted  against  their 
dynamic  pressure  counterparts  The  log-log  plot  in  Figure  17  exhibits  a  one-to-one  mapping, 
suggesting  that  the  dynamic  pressure  may  be  used  as  the  single  gain-scheduling  input  for 
the  fl  i  ght  envel  ope  of  consi derati  on. 
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"Q_op"  versus  "Dynamic  Pressure 
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Figurell:  Optimum  Q  versus  dynamic  pressure 


3.  Output  Feedback  Control  System 

3.1  Formulation 

In  the  earlier  stages  of  control  design,  we  assumed  that  all  the  states  were  accurately 
measureable.  This  assumption  does  not  hold  in  practice.  In  the  given  model,  6  out  of  the  24 
states  are  actually  measurable  (ie.,  IMU  outputs).  Therefore,  it  is  appropriate  to  include  a 
state  observer  design  as  shown  in  Figure  18. 

The  state  observer  estimates  all  24  states  including  the  six  measurable  ones  (ie.,  full  state 
observer)  by  adjusting  the  filter  gainKf  so  that  the  measurement  difference  (z  -  Jx )  remains 
small. 


Figure  18:  0  bserved-state  feedback  con  trol  system 
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In  order  to  derive  the  state  observer,  we  first  consider  the  pi  ant  and  measurement  equations 
x  =  Fx  +  Gu  (36) 

z  =  Jx .  (37) 

Suppose  we  construct  the  dynamic  equations  for  the  estimated  state  x  by  replicating  (36). 

x  =  Fx  +  Gu  (38) 

Now  define  the  estimation  error  as  ex  =  x-x,  and  derive  the  error  propagation  equation. 
Subtracting  (38)  from  (36),  we  arrive  at 

ex  =  x  -  x  =  Fx  -  Fx  =  Fex  .  (39) 

This  shows  that  when  the  matrix  F  is  asymptotically  stable,  the  error  ex converges  to  zero 
for  any  input  u,  which  means  that  x  eventually  converges  to  x  as  t  ->  oo .  However,  when 
F  is  not  stable,  ex  is  unbounded  and  x  may  grow  very  large.  To  avoid  this,  one  must  add  a 
correctional  term  to  (38)  as  shown  below. 

x  =  Fx  +  Gu  +  Kf(z- Jx)  (40) 

whereKf  isa24x6filter  gain  matrix. 

When  x  becomes  very  close  tox  ,  the  correctional  term  Kf(z- Jx)  plays  minimum  role. 
H  owever,  when  x  drifts  away  fromx ,  the  role  of  this  term  (driving  x  towards  x )  becomes 
more  dominant. 

To  see  how  this  can  be  done,  the  error  propagation  equation  is  given  by  subtracting  (40) 
from  (36), 

ex  =  x  -  x  =  (F  -  Kf  J)ex  (41) 

Now  ex  converges  to  zero  as  long  as  (F-KfJ)  is  asymptotically  stable.  Weknow  that  even 
if  F  is  unstable,  we  can  find  Kfthat  makes  (F-KfJ)  stable.  The  system  (40)  can  also  be 
expressed  as 


x  =  (F -Kf  J)x  + Gu +  Kfz  (42) 

This  state  observer  has  two  inputs  -  the  plant  input  uand  the  measurement  output  z  as 
shown  in  Figure  18. 

Any  choice  of  Kfin  (42)  for  which  (F-KfJ)is  asymptotically  stable  will  make  x  converge 
tox .  However,  in  general  the  output  z  and  the  process  dynamics  are  contaminated  by 
measurement  noise  and  disturbance  (process  noise)  respectively,  therefore  a  more  realistic 
model  representation  would  be 
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x  =  Fx  +  Gu  +  w  (43) 

z  =  Jx  +  v  (44) 

The  process  and  measurement  noises  w  and  v  are  assumed  Gaussian  and  satisfy 

E( w)  =  0  ,E(\)  -  0,is(wwr)  =  Qf ,  andisfvv^)  =  Rf ,  (45) 

Including  the  error  terms,  the  error  dynamics  in  (41)  can  be  rewritten  as 

ex  =Fx  +  w-(F-KfJ)x-Kf(Jx  +  v) 

=  (F-KfJ)ex  +w-Kfv.  (46) 


Because  of  w  and  v,  the  estimation  errors  may  not  converge  to  zero,  but  may  become  very 
small  if  the  rightKfis  selected.  This  motivates  the  so-called  Linear  Quadratic  Gaussian 
(LQG)  estimation  problem. 

The  obj  ecti  ve  i  s  to  fi  nd  K  f  that  d  ri  ves  the  fol  I  ow  i  ng  error  covari  ance  to  the  mi  ni  mu  m. 

Pf  =£(exe[),  (47) 

NotethatPf  is  symmetric  and  positive  semi-definite.  Solving  the  LQG  problem  leads  to  the 
differential  Riccati  equation  given  as 

Pf  =  FPf  +  PfFr  - Pf  JrR'‘jPf  +  Qf  (48) 

which  is  solved  forward  intimewith  theinitial  condition  of  Pf(r0) . 

If  the  observer  is  stable,  then  we  may  set  t0=- qo  ,  and  obtain  the  steady  state  solution  by 
solving  the  filtering  algebraic  Ricatti  equation  (FARE) 

0  =  FPf  +  PfFr  -  Pf  JrRf  JPf  +  Qf  (49) 

Using  the  analogy  of  this  equation  (FARE)  with  the  previous  CARE  in  the  LQR  design,  the 
analytic  method  of  [Naidu  2002]  or  matlab  function  care  can  be  used  again  to  generatePf . 

Pf  =care(  Fr,Jr,Qf,Rf) 
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And  the  optimal  filter  gain  Kf  is  given  as 

Kf=PfJrRf\  (50) 

3.2  Simulation 

Having  found  the  optimal  Kf,  the  state  observer  takes  the  measurements  and  the  control 
input,  and  estimates  the  states.  In  this  section,  we  repeat  the  first  two  examples  in  Section  2.4 
to  observe  the  effect  of  including  the  state  observer  on  the  performance.  The  simulation 
results  are  shown  in  Figure  19  and  Figure  20. 

The  errors  w  and  v  are  given  as  Gaussian  noises  [//  =  0, a  =  0.01]  and  their  covariance 
matrices  are  set  as 


Qf  —  1 24x24 

Rf  =10^1(5x6 

Figure  19  and  Figure  20  indicate  that  all  the  excited  states  are  estimated  well.  For  benign 
states  (ie.,  appearing  as  noises),  their  estimates  are  simply  the  noise-attenuated.  In  all  cases, 
no  estimate  breaks  aw  ay  from  the  true  states,  suggesting  that  F-KfJ  in  (46)  isstable.  Itwas 
also  observed  thatavery  highQf  /Rf  ratio  often  leads  to  instability. 


lime  (s)  time  (s) 


time  (s)  time  (s) 


- estimated  states 


■true  states 


time  (s)  time  (s)  time  (s)  time  (s) 


time  (s)  time  (s)  time  (s)  time  (s) 


Figure  19:  Plot  of  the  24  estimated  states  (blue)  versus  the  true  states(red)  for  the  high  altitude  low 
speed  run  with  0=10000 
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Figure20:  Plotofthe24  estimated  states  (blue)  and  the  true  states(red)  for  the  high  altitude  and  high 
speed  run  with  0=2. 


3.3  Effect  of  Adding  the  0  bserver  on  Closed  Loop  System 

The  step  responses  of  the  output  feedback  system  are  similar  to  those  of  the  state  feedback 
system.  It  appears  that  the  dynamical  behaviours  of  the  closed  loop  system  are  not  largely 
affected  by  inclusion  of  the  state  observer. 

Si  nee  the  dynamical  behaviour  is  solely  due  to  the  pole  locations,  we  can  ignore  the  input 
gain  terms  for  the  foil  owing  analysis.  We  revisit  the  process  equation  (36)  with  the  feedback 
input  u  =  -Kcx  and  also  the  estimation  error  differential  equation  (41). 


x  =  Fx  +  G(-Kcx) 

=  Fx-GKc(x-ex) 

=  (F  -  GKc)x  +  GKcex 

ex  =  Fx-Fx 
-  (F -Kf J)ex 


(51) 


(52) 


The  above  equations  can  be  combined  in  single  matrix  notation  as 
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F-GKC 

GKC 

X 

0 

F-KfJ 

_e,_ 

(53) 


The  above  equation  describes  the  dynamical  behaviour  of  the  output  feedback  control 
system.  The  characteristic  equation  becomes 


si  -  (F  -  GKC)  -GKC 

0  si  -  (F  -Kf  J) 


si  -  (F  -  GKc)||sI  -  (F  -  Kf  J)|  =  0  . 


(54) 


The  closed-loop  poles  of  the  combined  system  area  union  of  the  two  independent  groups  of 
poles  -  one  belonging  to  the  regulator  dynamics  F-GKcand  the  other  belonging  to  the 
state  observer  dynamicsF-KfJ  .  The  observer  does  not  shift  the  pole  locations  of  the 
regulator,  hence  making  it  possible  to  design  the  regulator  and  observer  separately  and 
bring  them  together  later.  This  property  is  known  as  the  Separation  Principle  [Anderson  and 
Moore  1989].  If  both  the  matrices  are  asymptotically  stable,  then  so  will  be  the  closed-loop 
system. 

From  the  steady  state  point  of  view,  inclusion  of  the  estimator  makes  no  difference 
(becausee  ->  0 ).  From  the  transient  point  of  view,  the  behaviour  of  the  combined  system 
might  be  different  from  that  of  the  true-state  feedback  system.  Ideally,  the  poles  of 
F-KfJ  should  be  placed  to  the  left  of  the  poles  of  F-GKcso  that  the  observer  response  is 
faster  than  the  regulator  response. 


The  current  setting  of  Rfand  Qf  results  in  the  observer  poles  placed  to  the  left  of  the 
regulator  poles  (dominant  poles  are  shown  below),  and  there  is  no  apparent  detrimental 
effects  from  the  observer  poles  on  the  overall  system  dynamics.  Further  experiments 
indicated  that  in  order  to  separate  the  observer  poles  further  from  the  regulator  poles,  the 
Qf  /Rf  ratio  needs  to  increase  further,  which  speeds  up  the  estimation  convergence.  As 
stated  before,  increasing  Qf  /Rf  excessively  leadsto  instability. 


Low  Speed 

Observer  poles 
-0.00018202 
-2.6  ±1.6e-014i 


Regulator  poles 
-5.8343^006 
-1.5082  ±1.5075i 


High  Speed 

Observer  poles 

-0.0010883 

-43.708 

-41.544  ±57.937i 


Regulator  poles 

-0.00032818 

-18.03 

-23.887  ±  31. 142i 
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4.  Conclusions 

We  present  an  optimal  flight  control  design  where  the  objective  is  to  find  a  feedback  gain 
that  achieves  a  good  balance  of  fast  response  and  small  control  effort.  We  adopt  the  generic 
24state  airframe  mod  el  derived  in  [Farhan  2007],  and  the  aerodynamic  coefficients  are  taken 
from  the  open  literature  [Garnell  and  East  1977]  where  the  missile  is  symmetric  and  has  a 
cruciform  shape. 

Initially,  effects  of  varying  the  weighting  matrix  Q  on  the  step  responses  were  investigated. 
It  was  concluded  that  increasing  Q  improves  the  response  speed  but  at  the  expense  of 
increased  control  burden.  The  closed  loop  responses  were  examined  at  numerous  operating 
points  along  the  flight  envelope,  and  the  corresponding  optimum  Q  values  are  noted.  When 
the  dynamic  pressure  becomes  high,  reducing  the  control  effort  becomes  the  primary 
objective  as  the  fast  response  is  readily  achievable.  It  looks  feasible  to  use  the  dynamic 
pressure  as  a  sol  e  gai  n-sched ul  i  ng  i  nput  for  the fl  i  ght  envel  ope  of  consi  derati  on . 

The  dominant  zeros  affect  the  initial  transient  of  the  step  response.  Tail -control led  missiles 
are  expected  to  exhibit  non-minimum  phase  responses.  However,  having  the  IMU  placed 
sufficiently  ahead  of  the  CG,  allows  the  positive  zero(s)  of  the  transfer  function 
(fly  out  /<3y  desired)  to  become  negative,  reversing  the  sign  of  the  initial  bump. 

The  design  and  simulation  of  an  LQG  controller  are  also  presented.  When  the  process  noise 
covariance  Qf  is  made  larger  than  the  measurement  noise  covarianceRf ,  the  estimator 
interprets  a  large  deviation  of  zfrom  z  as  an  indication  that  the  estimate  x  is  poor  hence 
needs  to  be  aggressively  corrected.  In  practice,  this  leads  to  large  Kf  and  fast  poles 
forF-KfJ  leading  to  fast  convergence.  Making  Qf  excessively  large,  however,  makes  the 
system  unstable.  Overall,  the  inclusion  of  the  state  observer  does  not  seem  to  deteriorate  the 
tracking  performance. 

The  LQR  based  controller  is  relatively  easy  to  implement  on  multivariable  systems,  and  the 
closed  loop  system  is  asymptotically  stable.  The  tuning  process  is  a  simple  adjustment  of 
Q/R  ratio  until  the  desired  balance  of  fast  response  and  small  control  effort  is  achieved. 
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